////////////////////////////////////////////////////////////////////////////////

*Install all necessary packages from SSC

////////////////////////////////////////////////////////////////////////////////
cap ssc install reghdfe
cap ssc install ftools
cap ssc install ppml
cap ssc install ppmlhdfe
cap ssc install mvdcmp
cap ssc install estout
cap ssc install mkest

////////////////////////////////////////////////////////////////////////////////

*Program to calculate marginal effects of year dummies in FE Poisson regressions

////////////////////////////////////////////////////////////////////////////////

cap program drop poisson_year_mfx

program define poisson_year_mfx
syntax [if] [in], save(string) [controls]
	qui {
	matrix define b=e(b)
	matrix define V=e(V)
*Generate predicted probabilities and SEs of BHPS responses
	predict linear_predict if insample==1, xb
	predict se if insample==1, stdp
	drop if linear_predict==.
	
	if "`controls'" != "" {
		if `"`if'"'!="" {
			local if_clause = `"`if'&"'
		}
		else if `"`if'"'=="" {
			local if_clause = "if "
		}
		foreach j in 1 2 3 {
			forval i=1991/2018 {
				local c = colnumb(matrix(b), "1.v`j'#`i'.year")
				scalar define b`j'_`i' = b[1,`c']
			}
		}
	
		foreach j in 1 2 3 {
			forval y=1991/2018 {
				local c = colnumb(matrix(V), "1.v`j'#`y'.year")
				local r = rownumb(matrix(V), "1.v`j'#`y'.year")
				scalar var`y'_v`j' = V[`r',`c']
			}
		}

*E[Y|X] by exponentiating the linear prediction-- *_0 is E[Y|X] in the base year
		gen exp_linear_predict=exp(linear_predict)
		gen exp_linear_predict_0=.
		gen linear_predict_0=.
		foreach j in 1 2 3 {
			forval i=1991/2017 {
				replace exp_linear_predict_0=exp(linear_predict-b`j'_`i') ///
				`if_clause't==`j'&year==`i'
				replace linear_predict_0=linear_predict-b`j'_`i' `if_clause't==`j'&year==`i'
			}
		}
		replace linear_predict_0=linear_predict `if_clause't==4	
		replace exp_linear_predict_0=exp_linear_predict `if_clause't==4
	
		gen numerator=exp_linear_predict
		gen numerator_0=exp_linear_predict_0
		bys person_year_id: egen denominator=total(numerator)
		bys person_year_id: egen denominator_0=total(numerator_0)
		gen probability_of_response=numerator/denominator
		gen probability_of_response_0=numerator_0/denominator_0

*Generate predicted BHPS responses
		gen x=probability_of_response*t
		bys person_year_id: egen expected_response=total(x)
		gen x_0=probability_of_response_0*t
		bys person_year_id: egen expected_response_0=total(x_0)
		
*Standard errors
		gen variance_ln_predict=se^2
		gen variance_ln_predict_0=.
		foreach i in 1 2 3 {
			forval y=1991/2017 {
				replace variance_ln_predict_0=var`y'_v`i'+se^2 `if_clause'year==`y'&t==`i'
			}
		}
		replace variance_ln_predict_0=se^2 `if_clause't==4
			

		forval y=1991/2018 {
			sum linear_predict `if_clause'year==`y'
			scalar define _mu_`y' = r(mean)
			scalar define N_`y' = r(N)
			sum linear_predict_0 `if_clause'year==`y'
			scalar define _mu_0_`y'=r(mean)
		}	

		gen variance_predict=.
		gen variance_predict_0=.
	
*Delta method for SEs
		gen g_prime_squared=(numerator*denominator-numerator^2)^2/denominator^4
		gen g_prime_squared_0= ///
		(numerator_0*denominator_0-numerator_0^2)^2/denominator_0^4
		gen variance_outcome=variance_ln_predict*g_prime_squared
		gen variance_outcome_0=variance_ln_predict_0*g_prime_squared_0
	
*Generate marginal effects of year dummies
		keep person_year_id probability_of_response* variance_outcome* bhps_security ///
		t year relative_xw
		drop if missing(probability_of_response_0)
		reshape wide probability_of_response probability_of_response_0 ///
		variance_outcome variance_outcome_0, i(person_year_id) j(t)
		collapse probability_of_response* variance_outcome* [pw=relative_xw], by(year)
		forval i=1/4 {
			forval y=1991/2018 {
				replace variance_outcome`i'=variance_outcome`i'/N_`y' if year==`y'
			}
		}	
	
*Marginal effect on probabilities--gives AMEs
		forval i=1/4 {
			gen dprob`i'dx_controls= ///
			probability_of_response`i'-probability_of_response_0`i'
		}
 
		rename variance_outcome* variance_outcome*_controls
		keep variance_outcome* dprob* year

		forval i=1/4 {
			rename variance_outcome`i'_controls variance`i'_controls
		}

		forval i=1/4 {
			rename variance_outcome_0`i'_controls variance`i'_controls_0
		}
	}
	
	else if "`controls'"=="" {
		gen exp_linear_predict=exp(linear_predict)
		gen numerator=exp_linear_predict
		bys person_year_id: egen denominator=total(numerator)
		gen probability_of_response=numerator/denominator

*Generate predicted BHPS responses
		gen x=probability_of_response*t
		bys person_year_id: egen expected_response=total(x)

*Standard errors
*Delta method
		gen variance_ln_predict=se^2
		gen g_prime_squared=(numerator*denominator-numerator^2)^2/denominator^4
		gen variance_outcome=variance_ln_predict*g_prime_squared


		forval y=1991/2018 {
			sum linear_predict if year==`y'
			scalar define N_`y' = r(N)
		}

*Generate marginal effects of year dummies
		keep person_year_id probability_of_response expected_response variance_outcome ///
		bhps_security t year relative_xw
		reshape wide probability_of_response variance_outcome, i(person_year_id) j(t)
		collapse probability* expected_response variance_outcome* [pw=relative_xw], by(year)
	

		forval i=1/4 {
			forval y=1991/2018 {
				replace variance_outcome`i'=variance_outcome`i'/N_`y' if year==`y'
			}
		}
	

*Marginal effect on probabilities
		forval i=1/4 { 
			gen probability_of_response`i'_2001=.
			replace probability_of_response`i'_2001=probability_of_response`i' ///
			if year==2001
			replace probability_of_response`i'_2001=0 if year!=2001
			egen probability_of_response`i'_2001_t=total(probability_of_response`i'_2001)
			replace probability_of_response`i'_2001=probability_of_response`i'_2001_t
			drop probability_of_response`i'_2001_t
			gen dprob`i'dx=probability_of_response`i'-probability_of_response`i'_2001
			gen variance_outcome`i'_2001=.
			replace variance_outcome`i'_2001=variance_outcome`i' if year==2001
			replace variance_outcome`i'_2001=0 if year!=2001
			egen variance_outcome`i'_2001_t=total(variance_outcome`i'_2001)
			replace variance_outcome`i'_2001=variance_outcome`i'_2001_t
			drop variance_outcome`i'_2001_t
			gen variance`i'_nocontrols=variance_outcome`i'+variance_outcome`i'_2001
		}

*Marginal effect on expected response	
		gen expected_response_2001=.
		replace expected_response_2001=expected_response if year==2001
		replace expected_response_2001=0 if year!=2001
		egen expected_response_2001_t=total(expected_response_2001)
		replace expected_response_2001=expected_response_2001_t
		drop expected_response_2001_t
		gen dexpected_responsedx=expected_response-expected_response_2001
	
		keep year variance1_nocontrols variance2_nocontrols variance3_nocontrols ///
		variance4_nocontrols dprob* variance_outcome*2001*
	
		forval i=1/4 {
			rename variance_outcome`i'_2001 variance`i'_2001
		}
	}
	}
	display "Year dummy mfx calculated. Saving `save'"
	save `save', replace
end

////////////////////////////////////////////////////////////////////////////////

*Program to calculate marginal effects of covariates in FE Poisson regressions

////////////////////////////////////////////////////////////////////////////////

cap program drop poisson_covariates_mfx

program define poisson_covariates_mfx
syntax varlist(min=1) [, continuous]
qui {
	foreach var in `varlist' {
		use "temp/poisson_with_predictions_controls.dta", clear
		gen exp_linear_predict=exp(linear_predict)
		gen exp_linear_predict_0=.
		gen linear_predict_0=.
		
		if "`continuous'"=="" {
			foreach j in 1 2 3 {
				local c = colnumb(matrix(b), "1.v`j'#1.`var'")
				scalar define b`j' = b[1,`c']
			}
	
			foreach j in 1 2 3 {
				local c = colnumb(matrix(V), "1.v`j'#1.`var'")
				local r = rownumb(matrix(V), "1.v`j'#1.`var'")
				scalar var_v`j' = V[`r',`c']
			}
			foreach j in 1 2 3 {
				replace exp_linear_predict_0=exp(linear_predict-b`j') if t==`j'&`var'==1
				replace linear_predict_0=linear_predict-b`j' if t==`j'&`var'==1
				replace exp_linear_predict_0=exp_linear_predict if t==`j'&`var'==0
				replace linear_predict_0=linear_predict if t==`j'&`var'==0
				replace exp_linear_predict=exp(linear_predict+b`j') if t==`j'&`var'==0
				replace linear_predict=linear_predict+b`j' if t==`j'&`var'==0
			}
	
			replace linear_predict_0=linear_predict if t==4&`var'==1	
			replace exp_linear_predict_0=exp_linear_predict if t==4&`var'==1
			replace linear_predict=linear_predict_0 if t==4&`var'==0	
			replace exp_linear_predict=exp_linear_predict_0 if t==4&`var'==0
		}
		else if "`continuous'"!="" {
			foreach j in 1 2 3 {
				local c = colnumb(matrix(b), "1.v`j'#c.`var'")
				scalar define b`j' = b[1,`c']
			}
	
			foreach j in 1 2 3 {
				local c = colnumb(matrix(V), "1.v`j'#c.`var'")
				local r = rownumb(matrix(V), "1.v`j'#c.`var'")
				scalar var_v`j' = V[`r',`c']
			}
			foreach j in 1 2 3 {
				replace exp_linear_predict_0=exp(linear_predict-b`j') if t==`j'
				replace linear_predict_0=linear_predict-b`j' if t==`j'
			}
			
			replace linear_predict_0=linear_predict if t==4	
			replace exp_linear_predict_0=exp_linear_predict if t==4
		}
	
		gen numerator=exp_linear_predict
		gen numerator_0=exp_linear_predict_0
		bys person_year_id: egen denominator=total(numerator)
		bys person_year_id: egen denominator_0=total(numerator_0)
		gen probability_of_response=numerator/denominator
		gen probability_of_response_0=numerator_0/denominator_0

*Generate predicted BHPS responses
		gen x=probability_of_response*t
		bys person_year_id: egen expected_response=total(x)
		gen x_0=probability_of_response_0*t
		bys person_year_id: egen expected_response_0=total(x_0)
		
*Standard errors
		gen variance_ln_predict=se^2
		gen variance_ln_predict_0=.
		foreach i in 1 2 3 {
			replace variance_ln_predict_0=var_v`i'+se^2 if t==`i'
		}
		replace variance_ln_predict_0=se^2 if t==4
			

		sum linear_predict
		scalar define _mu = r(mean)
		scalar define N = r(N)
		sum linear_predict_0
		scalar define _mu_0=r(mean)	

		gen variance_predict=.
		gen variance_predict_0=.
	
*Delta method
		gen g_prime_squared=(numerator*denominator-numerator^2)^2/denominator^4
		gen g_prime_squared_0= ///
		(numerator_0*denominator_0-numerator_0^2)^2/denominator_0^4
		gen variance_outcome=variance_ln_predict*g_prime_squared
		gen variance_outcome_0=variance_ln_predict_0*g_prime_squared_0
	
*Generate marginal effects of year dummies
		keep person_year_id probability_of_response* variance_outcome* bhps_security t relative_xw
		reshape wide probability_of_response probability_of_response_0 ///
		variance_outcome variance_outcome_0, i(person_year_id) j(t)
		collapse probability_of_response* variance_outcome* [pw=relative_xw]
		forval i=1/4 {
				replace variance_outcome`i'=variance_outcome`i'/N
		}	
	
*Marginal effect on probabilities
		forval i=1/4 {
			gen dprob`i'dx_controls= ///
			probability_of_response`i'-probability_of_response_0`i'
		} 
		rename variance_outcome* variance_outcome*_controls
		keep variance_outcome* dprob*
		forval i=1/4 {
			rename variance_outcome`i'_controls variance`i'_controls
		}
		forval i=1/4 {
			rename variance_outcome_0`i'_controls variance`i'_controls_0
		}

		gen covariate="`var'"
		gen ame= dprob1dx_controls+dprob2dx_controls
		gen se=sqrt(variance1_controls+variance2_controls+variance1_controls_0+variance2_controls_0)
		keep covariate ame se
	save "ames/`var'_new.dta", replace
}
}
end
